;;----------------------------------------------------------------------
;; scatter plot 
;;----------------------------------------------------------------------
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
  load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
  load "$KAILIBROOT/callvar.ncl" 
  load "$KAILIBROOT/callstr.ncl" 

  print("") 
  print("") 
  print("") 
  print("") 
  print("") 
  print("") 
;;----------------------------------------------------------------------
;; user setting 
;;----------------------------------------------------------------------

  colormap = "psgcap"
  icla = 148
  iclb = 106
  iclc = 159 ;;109
  icld = 178
  icle = 109

  LeftString   = ") Radon (mBq m~S~-3~N~ STP)"  
  RightString  = ""  
  XAxisString  = "Observation (mBq m~S~-3~N~ STP)~C~  "
  YAxisString  = "Nudged Simulation (mBq m~S~-3~N~ STP)"

  plotpath     = "~/fig/"
  plotpath     = "./" 
  plotname     = plotpath+"fig_scatter_region_grid"  
  plotform     = "eps"
  ploteps      = plotname+".eps"
  plotpdf      = plotname+".pdf"

  L_log = 0
  axisfont = 0.025 
  textfont = 0.022
  np = 12
  plot = new(np,graphic )

  markersize = 0.030 
  markertype = 1 

  abcd = (/"a","b","c","d","e","f","g","h","i","j","k","l","m","n"/) 

;;----------------------------------------------------------------------
;; read data 
;;----------------------------------------------------------------------
  ncfnama = "radon_mine_A0_1999-2003_ymm.nc" 
  ncfnamb = "radon_mine_A0_1999-2003_ymm.nc" 
  ncfnamc = "radon_mine_A2_1999-2003_ymm.nc" 
  ncfnamd = "radon_mine_A3_1999-2003_ymm.nc" 
  ncfname = "radon_obs_mm.nc" 
  ascnama = "radon_siteplot.dat" 
  ascnamb = "radon_sitename.dat" 
  ascnamc = "radon_sitebad.dat" 
  ascnamd = "radon_sitelat.dat" 
  ascname = "radon_sitelon.dat" 
  ascnamf = "radon_siteeur.dat" 
  ascnamg = "radon_sitechi.dat" 
  ascnamh = "radon_siteusa.dat" 
  ascnami = "radon_sitesam.dat" 

  fla = addfile(ncfnama,"r") 
  flb = addfile(ncfnamb,"r") 
  flc = addfile(ncfnamc,"r") 
  fld = addfile(ncfnamd,"r") 
  fle = addfile(ncfname,"r") 

  rma = fla->radon
  rmb = flb->radon
  rmc = flc->radon
  rmd = fld->radon
  rme = fle->radon

  kplot = asciiread(ascnama,-1,"integer") 
  sname = asciiread(ascnamb,-1,"string") 
  kbad  = asciiread(ascnamc,-1,"integer") 
  lat   = asciiread(ascnamd,-1,"float") 
  lon   = asciiread(ascname,-1,"float") 
  ieur  = asciiread(ascnamf,-1,"integer") 
  ichi  = asciiread(ascnamg,-1,"integer") 
  iusa  = asciiread(ascnamh,-1,"integer") 
  isam  = asciiread(ascnami,-1,"integer") 

  ns = dimsizes(kplot) 

  print("ns: "+ns) 

  do is = 0,ns-1
    if(ieur(is).eq.0) then 
       print("ieur: "+sname(is)) 
    end if 
  end do 

  do is = 0,ns-1
    if(ichi(is).eq.0) then 
       print("ichi: "+sname(is)) 
    end if 
  end do 

  do is = 0,ns-1
    if(isam(is).eq.0) then 
       print("isam: "+sname(is)) 
    end if 
  end do 

  do is = 0,ns-1
    if(iusa(is).eq.0) then 
       print("iusa: "+sname(is)) 
    end if 
  end do 


  knot  = kbad 
  kval  = kbad 


  rma!1 = "site" 
  rmb!1 = "site" 
  rmc!1 = "site" 
  rmd!1 = "site" 

;;----------------------------------------------------------------------
;; mask out bad sites 
;;----------------------------------------------------------------------
  do is = 0,ns-1 
     if(kbad(is).eq.0) then 
        rma(:,is) = -999  
        rmb(:,is) = -999  
        rmc(:,is) = -999  
        rmd(:,is) = -999  
        rme(:,is) = -999  
     end if  
  end do 

  do is = 0,ns-1
     if(all(ismissing(rme(:,is)))) then
        knot(is) = 0 
     end if
  end do
 
;;----------------------------------------------------------------------
;; annual average 
;;----------------------------------------------------------------------
  sma = dim_avg(rma(site|:,time|:))  
  smb = dim_avg(rmb(site|:,time|:))  
  smc = dim_avg(rmc(site|:,time|:))  
  smd = dim_avg(rmd(site|:,time|:))  
  sme = dim_avg(rme(site|:,time|:))  

  nt = 12 

;;----------------------------------------------------------------------
;; mask out sites only have annual mean 
;;----------------------------------------------------------------------

  do is = 0,ns-1

     xrmaa = rma(0,is) + rma(1,is) + rma(11,is)  
     xrmab = rma(2,is) + rma(3,is) + rma(4,is)  
     xrmac = rma(5,is) + rma(6,is) + rma(7,is)  
     xrmad = rma(8,is) + rma(9,is) + rma(10,is)  

     xrmba = rmb(0,is) + rmb(1,is) + rmb(11,is)  
     xrmbb = rmb(2,is) + rmb(3,is) + rmb(4,is)  
     xrmbc = rmb(5,is) + rmb(6,is) + rmb(7,is)  
     xrmbd = rmb(8,is) + rmb(9,is) + rmb(10,is)  

     xrmca = rmc(0,is) + rmc(1,is) + rmc(11,is)  
     xrmcb = rmc(2,is) + rmc(3,is) + rmc(4,is)  
     xrmcc = rmc(5,is) + rmc(6,is) + rmc(7,is)  
     xrmcd = rmc(8,is) + rmc(9,is) + rmc(10,is)  

     xrmda = rmd(0,is) + rmd(1,is) + rmd(11,is)  
     xrmdb = rmd(2,is) + rmd(3,is) + rmd(4,is)  
     xrmdc = rmd(5,is) + rmd(6,is) + rmd(7,is)  
     xrmdd = rmd(8,is) + rmd(9,is) + rmd(10,is)  

  ;if(is.eq.0) then    ;; beijing
  ;   rme(:,is) =  -999 
  ;   rme(0,is) =  9078. 
  ;   rme(3,is) =  6248. 
  ;   rme(6,is) =  6847. 
  ;   rme(9,is) =  9945. 
  ;end if  

  if(is.eq.93 .or. \ 
     is.eq.94 .or. \ 
     is.eq.95 .or. \ 
     is.eq.96 .or. \ 
     is.eq.97 .or. \ 
     is.eq.98 .or. \ 
     is.eq.100 ) then 
     rma(:,is) =  -999
     rmb(:,is) =  -999
     rmc(:,is) =  -999
     rmd(:,is) =  -999
     rma(0,is) = xrmaa / 3.  
     rmb(0,is) = xrmba / 3. 
     rmc(0,is) = xrmca / 3. 
     rmd(0,is) = xrmda / 3. 
     rma(3,is) = xrmab / 3. 
     rmb(3,is) = xrmbb / 3. 
     rmc(3,is) = xrmcb / 3. 
     rmd(3,is) = xrmdb / 3. 
     rma(6,is) = xrmac / 3. 
     rmb(6,is) = xrmbc / 3. 
     rmc(6,is) = xrmcc / 3. 
     rmd(6,is) = xrmdc / 3. 
     rma(9,is) = xrmad / 3. 
     rmb(9,is) = xrmbd / 3. 
     rmc(9,is) = xrmcd / 3. 
     rmd(9,is) = xrmdd / 3. 
  end if 

  if(is.eq.93) then   ;; huhehaote
     rme(:,is) =  -999 
     rme(0,is) = 10600.
     rme(3,is) =  6900.
     rme(6,is) =-10000.
     rme(9,is) =  7340. 
  end if   
  if(is.eq.94) then    ;; changchun
     rme(:,is) =  -999 
    ;rme(0,is) = 10500.
     rme(0,is) = 13000.
     rme(3,is) =  8800.
     rme(6,is) = 11380.
     rme(9,is) =  9390. 
  end if  
  if(is.eq.95) then   ;; nanjing  
     rme(:,is) =  -999
     rme(0,is) =  6410.
     rme(3,is) =  6420.
     rme(6,is) =  6430.
     rme(9,is) =  6440.
  end if
  if(is.eq.96) then   ;; xi_an
     rme(:,is) =  -999 
     rme(0,is) = 17400. 
     rme(3,is) =  9600. 
     rme(6,is) =  9200. 
     rme(9,is) = 10000
  end if 
  if(is.eq.97) then   ;; shanghai 
     rme(:,is) =  -999
     rme(0,is) =  4810.
     rme(3,is) =  4820.
     rme(6,is) =  4830.
     rme(9,is) =  4840.
  end if
  if(is.eq.98) then    ;; wuhan
     rme(:,is) =  -999 
     rme(0,is) = 23920.
     rme(3,is) = 13770.
     rme(6,is) =  8500.
     rme(9,is) = 12670. 
  end if  
  if(is.eq.100) then   ;; guiyang
     rme(:,is) =  -999
     rme(0,is) = 12210.
     rme(3,is) = 12220.
     rme(6,is) = 12230.
     rme(9,is) = 12240.
  end if

  if(is.eq.30) then   ;; rio
     rme(1,is) = 1450.
     rme(8,is) = 1560.
     rme(9,is) = 1460.
     rme(11,is)= 1390.
     print(sname(is)+" "+rme(:,is)) 
  end if

  if(is.eq.29) then   ;; cha
     rme(1,is) = 1450.
     rme(3,is) = 1490.
     rme(4,is) = 1470.
     rme(5,is) = 1600.
     rme(7,is) = 1800.
     rme(9,is) = 1700.
     rme(10,is) = 1680.
     print(sname(is)+" "+rme(:,is))
  end if

  end do 


  print("--------- valid ann sites -----------") 
  do is = 0,ns-1 
    ;print("-") 
    ;print("obs: "+sname(is)+"  "+rme(:,is)+"   "+rma(:,is)+"   "+rmd(:,is)) 
    ;print("-") 
    if(.not.ismissing(sme(is)) .and. .not. ismissing(rme(0,is)) ) then 
       if(abs(sme(is)-rme(0,is)).le.1.) then 
          rme(:,is) = -999 
          rma(:,is) = -999 
          rmb(:,is) = -999 
          rmc(:,is) = -999 
          rmd(:,is) = -999 
          knot(is)  = 0
          if(kbad(is).eq.1) then  
             print(strnice(sname(is),20)+"") 
          end if 
          kbad(is) = 0 
       end if 
    end if 
  end do 
  print("") 
  print("") 

  print("--------- valid sites -----------")
  do is = 0,ns-1 
    if(knot(is).ne.0) then 
       print(strnice(sname(is),20)+"       "+sprintf("%9.3f",sme(is))+"    "+sprintf("%9.3f",smd(is)) ) 
    else
       kbad(is) = 0 
    end if 
  end do 
  print("") 
  print("") 
  print("number of valid sites: "+sprinti("%6.0i",sum(kbad)) ) 
  print("") 
  print("") 
  print("------------------------------------------------------------------") 

;;----------------------------------------------------------------------
;; recalculate annual mean 
;;----------------------------------------------------------------------
  sma = dim_avg(rma(site|:,time|:))  
  smb = dim_avg(rmb(site|:,time|:))  
  smc = dim_avg(rmc(site|:,time|:))  
  smd = dim_avg(rmd(site|:,time|:))  
  sme = dim_avg(rme(site|:,time|:))  
 
;;----------------------------------------------------------------------
;; arrage variables  
;;----------------------------------------------------------------------

  rmea = rme 
  rmea1= rme 
  rmea2= rme 
  rmeb = rme 
  rmec = rme 
  rmed = rme 

  rmaa = rma 
  rmaa1= rma 
  rmaa2= rma 
  rmab = rma 
  rmac = rma 
  rmad = rma 

  rmca = rmc 
  rmca1= rmc 
  rmca2= rmc 
  rmcb = rmc 
  rmcc = rmc 
  rmcd = rmc 

  rmda = rmd 
  rmda1= rmd 
  rmda2= rmd 
  rmdb = rmd
  rmdc = rmd 
  rmdd = rmd 

  do is = 0,ns-1 
  if(ieur(is).ne.0) then 
     rmea(:,is) = -999 
     rmaa(:,is) = -999 
     rmca(:,is) = -999 
     rmda(:,is) = -999 
     rmea1(:,is) = -999
     rmaa1(:,is) = -999
     rmca1(:,is) = -999
     rmda1(:,is) = -999
     rmea2(:,is) = -999
     rmaa2(:,is) = -999
     rmca2(:,is) = -999
     rmda2(:,is) = -999
  end if

  if(ieur(is).eq.0) then 
  if(lat(is).gt.60.) then
     rmea1(:,is) = -999 
     rmaa1(:,is) = -999 
     rmca1(:,is) = -999 
     rmda1(:,is) = -999 
  else
     rmea2(:,is) = -999
     rmaa2(:,is) = -999
     rmca2(:,is) = -999
     rmda2(:,is) = -999
  end if
  end if

  print(sname(is)+"   "+max(rme(:,is))) 
  end do 

  do is = 0,ns-1
  if(ichi(is).ne.0) then
     rmeb(:,is) = -999
     rmab(:,is) = -999
     rmcb(:,is) = -999
     rmdb(:,is) = -999
  end if
  print(sname(is)+"   "+max(rme(:,is)))
  end do

  do is = 0,ns-1
  if(iusa(is).ne.0) then
     rmec(:,is) = -999
     rmac(:,is) = -999
     rmcc(:,is) = -999
     rmdc(:,is) = -999
  end if
  print(sname(is)+"   "+max(rme(:,is)))
  end do

  do is = 0,ns-1
  if(isam(is).ne.0) then
     rmed(:,is) = -999
     rmad(:,is) = -999
     rmcd(:,is) = -999
     rmdd(:,is) = -999
  end if
  print(sname(is)+"   "+max(rme(:,is)))
  end do

  x1  = ndtooned(rme(:,:)) 
  y1  = ndtooned(rma(:,:)) 
  y2  = ndtooned(rmc(:,:)) 
  y3  = ndtooned(rmd(:,:)) ;; a shift !!!  

  x1a  = ndtooned(rmea(:,:)) 
  y1a  = ndtooned(rmaa(:,:)) 
  y2a  = ndtooned(rmca(:,:)) 
  y3a  = ndtooned(rmda(:,:)) ;; a shift !!!  

  x1a1  = ndtooned(rmea1(:,:))
  y1a1  = ndtooned(rmaa1(:,:))
  y2a1  = ndtooned(rmca1(:,:))
  y3a1  = ndtooned(rmda1(:,:)) ;; a shift !!!

  x1a2  = ndtooned(rmea2(:,:))
  y1a2  = ndtooned(rmaa2(:,:))
  y2a2  = ndtooned(rmca2(:,:))
  y3a2  = ndtooned(rmda2(:,:)) ;; a shift !!!

  x1b  = ndtooned(rmeb(:,:)) 
  y1b  = ndtooned(rmab(:,:)) 
  y2b  = ndtooned(rmcb(:,:)) 
  y3b  = ndtooned(rmdb(:,:)) ;; a shift !!!  

  x1c  = ndtooned(rmec(:,:)) 
  y1c  = ndtooned(rmac(:,:)) 
  y2c  = ndtooned(rmcc(:,:)) 
  y3c  = ndtooned(rmdc(:,:)) ;; a shift !!!  

  x1d  = ndtooned(rmed(:,:)) 
  y1d  = ndtooned(rmad(:,:)) 
  y2d  = ndtooned(rmcd(:,:)) 
  y3d  = ndtooned(rmdd(:,:)) ;; a shift !!!  


;;----------------------------------------------------------------------
;; mask negative values again 
;;----------------------------------------------------------------------
  x1  = var_maskneg(x1,0) 
  y1  = var_maskneg(y1,0) 
  y2  = var_maskneg(y2,0) 
  y3  = var_maskneg(y3,0) 

  x1a  = var_maskneg(x1a,0) 
  y1a  = var_maskneg(y1a,0) 
  y2a  = var_maskneg(y2a,0) 
  y3a  = var_maskneg(y3a,0) 

  x1a1  = var_maskneg(x1a1,0)
  y1a1  = var_maskneg(y1a1,0)
  y2a1  = var_maskneg(y2a1,0)
  y3a1  = var_maskneg(y3a1,0)

  x1a2  = var_maskneg(x1a2,0)
  y1a2  = var_maskneg(y1a2,0)
  y2a2  = var_maskneg(y2a2,0)
  y3a2  = var_maskneg(y3a2,0)

  x1b  = var_maskneg(x1b,0) 
  y1b  = var_maskneg(y1b,0) 
  y2b  = var_maskneg(y2b,0) 
  y3b  = var_maskneg(y3b,0) 

  x1c  = var_maskneg(x1c,0) 
  y1c  = var_maskneg(y1c,0) 
  y2c  = var_maskneg(y2c,0) 
  y3c  = var_maskneg(y3c,0) 

  x1d  = var_maskneg(x1d,0) 
  y1d  = var_maskneg(y1d,0) 
  y2d  = var_maskneg(y2d,0) 
  y3d  = var_maskneg(y3d,0) 

  ;print("x1,y1,y2: "+x1+"   "+y1+"   "+y2) 

  fy1 = var_p2(x1,y1,0) 
  fy2 = var_p2(x1,y2,0) 
  fy3 = var_p2(x1,y3,0) 
  
  fy1a = var_p2(x1a,y1a,0) 
  fy2a = var_p2(x1a,y2a,0) 
  fy3a = var_p2(x1a,y3a,0) 
  
  fy1a1 = var_p2(x1a1,y1a1,0)
  fy2a1 = var_p2(x1a1,y2a1,0)
  fy3a1 = var_p2(x1a1,y3a1,0)

  fy1a2 = var_p2(x1a2,y1a2,0)
  fy2a2 = var_p2(x1a2,y2a2,0)
  fy3a2 = var_p2(x1a2,y3a2,0)
 
  fy1b = var_p2(x1b,y1b,0) 
  fy2b = var_p2(x1b,y2b,0) 
  fy3b = var_p2(x1b,y3b,0) 
  
  fy1c = var_p2(x1c,y1c,0) 
  fy2c = var_p2(x1c,y2c,0) 
  fy3c = var_p2(x1c,y3c,0) 
  
  fy1d = var_p2(x1d,y1d,0) 
  fy2d = var_p2(x1d,y2d,0) 
  fy3d = var_p2(x1d,y3d,0) 
  
  print(""+fy1+" "+fy2+" "+fy3)
 
;;----------------------------------------------------------------------
;; plotvar  
;;----------------------------------------------------------------------
  yy = new((/3,dimsizes(y1)/),"float") 

  yy(0,:) = y1(:) 
  yy(1,:) = y2(:) 
  yy(2,:) = y3(:) 

  ya = new((/3,dimsizes(y1)/),"float") 
  ya(0,:) = y1a(:) 
  ya(1,:) = y2a(:) 
  ya(2,:) = y3a(:) 

  ya1 = new((/3,dimsizes(y1)/),"float")
  ya1(0,:) = y1a1(:)
  ya1(1,:) = y2a1(:)
  ya1(2,:) = y3a1(:)


  ya2 = new((/3,dimsizes(y1)/),"float")
  ya2(0,:) = y1a2(:)
  ya2(1,:) = y2a2(:)
  ya2(2,:) = y3a2(:)

  yb = new((/3,dimsizes(y1)/),"float") 
  yb(0,:) = y1b(:) 
  yb(1,:) = y2b(:) 
  yb(2,:) = y3b(:) 

  yc = new((/3,dimsizes(y1)/),"float") 
  yc(0,:) = y1c(:) 
  yc(1,:) = y2c(:) 
  yc(2,:) = y3c(:) 

  yd = new((/3,dimsizes(y1)/),"float") 
  yd(0,:) = y1d(:) 
  yd(1,:) = y2d(:) 
  yd(2,:) = y3d(:) 


;;----------------------------------------------------------------------
;; number of non-missing samples  
;;----------------------------------------------------------------------
  nxnmisa = where(ismissing(x1a),0,1) 
  nxnmisb = where(ismissing(x1b),0,1) 
  nxnmisc = where(ismissing(x1c),0,1) 
  nxnmisd = where(ismissing(x1d),0,1) 

  nnmisa  = dim_sum(nxnmisa)
  nnmisb  = dim_sum(nxnmisb)  
  nnmisc  = dim_sum(nxnmisc)  
  nnmisd  = dim_sum(nxnmisd)  

;;----------------------------------------------------------------------
;; correlation 
;;----------------------------------------------------------------------
  fr1 = esccr(x1,y1,0)
  fr2 = esccr(x1,y2,0)
  fr3 = esccr(x1,y3,0)
 
  fr1a = esccr(x1a,y1a,0)
  fr2a = esccr(x1a,y2a,0)
  fr3a = esccr(x1a,y3a,0)

  fr1a1 = esccr(x1a1,y1a1,0)
  fr2a1 = esccr(x1a1,y2a1,0)
  fr3a1 = esccr(x1a1,y3a1,0)

  fr1a2 = esccr(x1a2,y1a2,0)
  fr2a2 = esccr(x1a2,y2a2,0)
  fr3a2 = esccr(x1a2,y3a2,0)
 
  fr1b = esccr(x1b,y1b,0)
  fr2b = esccr(x1b,y2b,0)
  fr3b = esccr(x1b,y3b,0)
 
  fr1c = esccr(x1c,y1c,0)
  fr2c = esccr(x1c,y2c,0)
  fr3c = esccr(x1c,y3c,0)
 
  fr1d = esccr(x1d,y1d,0)
  fr2d = esccr(x1d,y2d,0)
  fr3d = esccr(x1d,y3d,0)
 
  print("") 
  print(""+fr1+" "+fr2+" "+fr3)


  cr = new(np,"float") 
  p2 = new(np,"float") 
  nn = new(np,"integer") 
  cr = 0. 
  p2 = 0. 
  nn = 0

  cr(0) = fr1a
  cr(1) = fr2a
  cr(2) = fr3a
  cr(3) = fr1b 
  cr(4) = fr2b 
  cr(5) = fr3b 
  cr(6) = fr1c 
  cr(7) = fr2c 
  cr(8) = fr3c 
  cr(9) = fr1d 
  cr(10)= fr2d 
  cr(11)= fr3d 

  p2(0) = fy1a
  p2(1) = fy2a
  p2(2) = fy3a
  p2(3) = fy1b
  p2(4) = fy2b
  p2(5) = fy3b
  p2(6) = fy1c
  p2(7) = fy2c
  p2(8) = fy3c
  p2(9) = fy1d
  p2(10)= fy2d
  p2(11)= fy3d

;;----------------------------------------------------------------------
;; open wks and plot 
;;----------------------------------------------------------------------
 

  wks = gsn_open_wks(plotform,plotname)
  gsn_define_colormap(wks,colormap)

;;----------------------------------------------------------------------
;; plot setting 
;;----------------------------------------------------------------------
  res          = True
  res@gsnDraw  = False
  res@gsnFrame = False

  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = YAxisString
  res@tiXAxisFontHeightF     = axisfont 
  res@tiYAxisFontHeightF     = axisfont 
  res@tmYLLabelFontHeightF   = axisfont
  res@tmXBLabelFontHeightF   = axisfont

  res@gsnLeftString              = LeftString
  res@gsnRightString             = RightString 
  res@gsnLeftStringFontHeightF   = axisfont 
  res@gsnRightStringFontHeightF  = axisfont 
 
 
  ;;res@mpShapeMode  = "FreeAspect"
  ;;res@vpWidthF      = 0.6
  ;;res@vpHeightF     = 0.6

  
  if(L_log.eq.0) 
    res@trXLog     = True
    res@trYLog     = True
    res@trXMinF    = 50
    res@trXMaxF    = 50000
    res@trYMinF    = 50
    res@trYMaxF    = 50000
  else
    res@trXLog     = False
    res@trYLog     = False
    res@trXMinF    = 0.
    res@trXMaxF    = 5000
    res@trYMinF    = 0.
    res@trYMaxF    = 5000
  end if
  
  res@trXReverse     = False

;;----------------------------------------------------------------------
;; marker setting 
;;----------------------------------------------------------------------
  res@xyMarkLineModes   = "Markers"
  res@xyMarkers         = markertype
  res@xyMarkerSizeF     = markersize
  res@xyMarkerColor     = icla 

;;----------------------------------------------------------------------
;; raw plot 
;;----------------------------------------------------------------------

  ip = 0 

  res@tiXAxisString       = "" 
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = YAxisString
  res@gsnLeftString       = "" ;;abcd(ip) + ") WCRP1995" 
  res@gsnRightString      = "" ;;"Europe" 

  res@xyMarkerColor     = icla
  res@xyMarkers         = markertype
  res@xyMarkerSizeF     = markersize
  res@xyMarkerThicknessF= 1.0 
  plot(ip) = gsn_csm_xy(wks,x1a1,ya1(0,:),res)

  res@tiXAxisString       = "" 
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = YAxisString
  res@gsnLeftString       = abcd(ip) + ") WCRP1995" 
  res@gsnRightString      = "Europe" 

  res@xyMarkerColor     = icle
  res@xyMarkers         = 2
  res@xyMarkerSizeF     = 0.01
  res@xyMarkerThicknessF= 0.5 

  plota    = gsn_csm_xy(wks,x1a2,ya2(0,:),res)

  overlay(plot(ip),plota) 


  ip = ip + 1 
  res@tiXAxisString       = "" 
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = "" ;;YAxisString
  res@gsnLeftString       = "" ;;abcd(ip) + ") WCRP1995" 
  res@gsnRightString      = "" ;;"Europe" 

  res@xyMarkerColor     = icla
  res@xyMarkers         = markertype
  res@xyMarkerSizeF     = markersize
  res@xyMarkerThicknessF= 1.0 
  plot(ip) = gsn_csm_xy(wks,x1a1,ya1(1,:),res)

  res@tiXAxisString       = "" 
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = "" 
  res@gsnLeftString       = abcd(ip) + ") SW1998 scaled"
  res@gsnRightString      = "Europe" 

  res@xyMarkerColor     = icle
  res@xyMarkers         = 2
  res@xyMarkerSizeF     = 0.01
  res@xyMarkerThicknessF= 0.5 
  plota    = gsn_csm_xy(wks,x1a2,ya2(1,:),res)

  overlay(plot(ip),plota) 

  ip = ip + 1 
  res@tiXAxisString       = "" 
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = "" ;;YAxisString
  res@gsnLeftString       = "" ;;abcd(ip) + ") WCRP1995" 
  res@gsnRightString      = "" ;;"Europe" 

  res@xyMarkerColor     = icla
  res@xyMarkers         = markertype
  res@xyMarkerSizeF     = markersize
  res@xyMarkerThicknessF= 1.0 
  plot(ip) = gsn_csm_xy(wks,x1a1,ya1(2,:),res)

  res@tiXAxisString       = "" 
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = "" 
  res@gsnLeftString       = abcd(ip) + ") Merged" 
  res@gsnRightString      = "Europe" 

  res@xyMarkerColor     = icle
  res@xyMarkers         = 2
  res@xyMarkerSizeF     = 0.01
  res@xyMarkerThicknessF= 0.5 
  plota    = gsn_csm_xy(wks,x1a2,ya2(2,:),res)

  overlay(plot(ip),plota) 

;;; 

  res@xyMarkers         = markertype
  res@xyMarkerSizeF     = markersize

  res@xyMarkerColor     = iclb 

  ip = ip + 1 
  res@tiXAxisString       = "" 
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = YAxisString
  res@gsnLeftString       = abcd(ip) + ") WCRP1995" 
  res@gsnRightString      = "China" 
  plot(ip) = gsn_csm_xy(wks,x1b,yb(0,:),res)

  ip = ip + 1 
  res@tiXAxisString       = "" 
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = ""
  res@gsnLeftString       = abcd(ip) + ") SW1998 scaled"
  res@gsnRightString      = "China" 
  plot(ip) = gsn_csm_xy(wks,x1b,yb(1,:),res)

  ip = ip + 1 
  res@tiXAxisString       = "" 
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = ""
  res@gsnLeftString       = abcd(ip) + ") Merged" 
  res@gsnRightString      = "China" 
  plot(ip) = gsn_csm_xy(wks,x1b,yb(2,:),res)

;;; 
  res@xyMarkerColor     = iclc 

  ip = ip+1
  res@tiXAxisString       = "" ;;XAxisString
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = YAxisString
  res@gsnLeftString       = abcd(ip) + ") WCRP1995" 
  res@gsnRightString      = "USA" 
  plot(ip) = gsn_csm_xy(wks,x1c,yc(0,:),res)

  ip = ip+1
  res@tiXAxisString       = "" ;XAxisString
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = ""
  res@gsnLeftString       = abcd(ip) + ") SW1998 scaled"
  res@gsnRightString      = "USA" 
  plot(ip) = gsn_csm_xy(wks,x1c,yc(1,:),res)

  ip = ip+1
  res@tiXAxisString       = "" ;;XAxisString
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = ""
  res@gsnLeftString       = abcd(ip) + ") Merged" 
  res@gsnRightString      = "USA" 
  plot(ip) = gsn_csm_xy(wks,x1c,yc(2,:),res)


;;;
  res@xyMarkerColor     = icld 

  ip = ip+1
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = YAxisString
  res@gsnLeftString       = abcd(ip) + ") WCRP1995"
  res@gsnRightString      = "South America"
  plot(ip) = gsn_csm_xy(wks,x1d,yd(0,:),res)

  ip = ip+1
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = ""
  res@gsnLeftString       = abcd(ip) + ") SW1998 scaled"
  res@gsnRightString      = "South America"
  plot(ip) = gsn_csm_xy(wks,x1d,yd(1,:),res)

  ip = ip+1
  res@tiXAxisString       = XAxisString
  res@tiYAxisString       = ""
  res@gsnLeftString       = abcd(ip) + ") Merged"
  res@gsnRightString      = "South America"
  plot(ip) = gsn_csm_xy(wks,x1d,yd(2,:),res)

;;----------------------------------------------------------------------
;; print P2 and CR values on plot 
;;----------------------------------------------------------------------
  txres               = True
  txres@txFontHeightF = textfont 
  txres@txJust        = "CenterLeft" 

  textc = new(np,graphic)  
  texta = new(np,graphic)  
  textb = new(np,graphic)  

  nn(0:2) = nnmisa 
  nn(3:5) = nnmisb 
  nn(6:8) = nnmisc 
  nn(9:11)= nnmisd 

  do ip = 0,np-1
     ;if(ip.eq.0.or.ip.eq.3.or.ip.eq.6) then 
     textc(ip) = gsn_add_text(wks,plot(ip),"  N = "+sprinti("%3.0i",nn(ip)),    110.,10000.,txres) 
     texta(ip) = gsn_add_text(wks,plot(ip), "P2 = "+sprintf("%3.2f",p2(ip))+"%",3000.,250.,txres)
     textb(ip) = gsn_add_text(wks,plot(ip),"  R = "+sprintf("%3.2f",cr(ip)),    3000.,150.,txres)
  end do 

;;----------------------------------------------------------------------
;; polygon 
;;----------------------------------------------------------------------
  xp    = new(6,"float")
  yp    = new(6,"float")

  xp(0) = 1.0 * res@trXMinF  
  yp(0) = 1.0 * res@trYMinF  
  xp(1) = 2.0 * res@trXMinF 
  yp(1) = 1.0 * res@trYMinF  
  xp(2) = 1.0 * res@trXMaxF 
  yp(2) = 0.5 * res@trYMaxF 
  xp(3) = 1.0 * res@trXMaxF 
  yp(3) = 1.0 * res@trYMaxF 
  xp(4) = 0.5 * res@trXMaxF  
  yp(4) = 1.0 * res@trYMaxF  
  xp(5) = 1.0 * res@trXMinF  
  yp(5) = 2.0 * res@trYMinF  

 ;gsres                   = True                        ; poly res
 ;gsres@tfPolyDrawOrder   = "Predraw"                   ; draw this first
 ;gsres@gsFillColor       = "grey"                      ; color chosen
 ;duma = gsn_add_polygon(wks,plot0,xp,yp,gsres)

;;----------------------------------------------------------------------
;; polyline 
;;----------------------------------------------------------------------
  resl = True
  resl@gsLineColor = "black"
  resl@gsLineLabelString= ""

  xpts = (/1.,1./)
  ypts = (/1.,1./)

  ypts(0) = res@trYMinF
  ypts(1) = res@trYMaxF
  xpts(0) = res@trXMinF
  xpts(1) = res@trXMaxF

  resl@gsLineThicknessF = 0.1
  resl@gsLineDashPattern = 0

  dunaa=gsn_add_polyline(wks,plot(0),xpts(0:1),ypts(0:1),resl)
  dunab=gsn_add_polyline(wks,plot(1),xpts(0:1),ypts(0:1),resl)
  dunac=gsn_add_polyline(wks,plot(2),xpts(0:1),ypts(0:1),resl)
  dunad=gsn_add_polyline(wks,plot(3),xpts(0:1),ypts(0:1),resl)
  dunae=gsn_add_polyline(wks,plot(4),xpts(0:1),ypts(0:1),resl)
  dunaf=gsn_add_polyline(wks,plot(5),xpts(0:1),ypts(0:1),resl)
  dunag=gsn_add_polyline(wks,plot(6),xpts(0:1),ypts(0:1),resl)
  dunah=gsn_add_polyline(wks,plot(7),xpts(0:1),ypts(0:1),resl)
  dunai=gsn_add_polyline(wks,plot(8),xpts(0:1),ypts(0:1),resl)
  dunaj=gsn_add_polyline(wks,plot(9),xpts(0:1),ypts(0:1),resl)
  dunak=gsn_add_polyline(wks,plot(10),xpts(0:1),ypts(0:1),resl)
  dunal=gsn_add_polyline(wks,plot(11),xpts(0:1),ypts(0:1),resl)

  resl@gsLineThicknessF = 1.0
  resl@gsLineDashPattern = 11

  ypts(0) = res@trYMinF
  ypts(1) = 0.5 * res@trYMaxF
  xpts(0) = 2. * res@trXMinF
  xpts(1) = res@trXMaxF

  dunba=gsn_add_polyline(wks,plot(0),xpts(0:1),ypts(0:1),resl)
  dunbb=gsn_add_polyline(wks,plot(1),xpts(0:1),ypts(0:1),resl)
  dunbc=gsn_add_polyline(wks,plot(2),xpts(0:1),ypts(0:1),resl)
  dunbd=gsn_add_polyline(wks,plot(3),xpts(0:1),ypts(0:1),resl)
  dunbe=gsn_add_polyline(wks,plot(4),xpts(0:1),ypts(0:1),resl)
  dunbf=gsn_add_polyline(wks,plot(5),xpts(0:1),ypts(0:1),resl)
  dunbg=gsn_add_polyline(wks,plot(6),xpts(0:1),ypts(0:1),resl)
  dunbh=gsn_add_polyline(wks,plot(7),xpts(0:1),ypts(0:1),resl)
  dunbi=gsn_add_polyline(wks,plot(8),xpts(0:1),ypts(0:1),resl)
  dunbj=gsn_add_polyline(wks,plot(9),xpts(0:1),ypts(0:1),resl)
  dunbk=gsn_add_polyline(wks,plot(10),xpts(0:1),ypts(0:1),resl)
  dunbl=gsn_add_polyline(wks,plot(11),xpts(0:1),ypts(0:1),resl)

  ypts(0) = 2. * res@trYMinF
  ypts(1) = res@trYMaxF
  xpts(0) = res@trXMinF
  xpts(1) = 0.5 * res@trXMaxF

  dunca=gsn_add_polyline(wks,plot(0),xpts(0:1),ypts(0:1),resl)
  duncb=gsn_add_polyline(wks,plot(1),xpts(0:1),ypts(0:1),resl)
  duncc=gsn_add_polyline(wks,plot(2),xpts(0:1),ypts(0:1),resl)
  duncd=gsn_add_polyline(wks,plot(3),xpts(0:1),ypts(0:1),resl)
  dunce=gsn_add_polyline(wks,plot(4),xpts(0:1),ypts(0:1),resl)
  duncf=gsn_add_polyline(wks,plot(5),xpts(0:1),ypts(0:1),resl)
  duncg=gsn_add_polyline(wks,plot(6),xpts(0:1),ypts(0:1),resl)
  dunch=gsn_add_polyline(wks,plot(7),xpts(0:1),ypts(0:1),resl)
  dunci=gsn_add_polyline(wks,plot(8),xpts(0:1),ypts(0:1),resl)
  duncj=gsn_add_polyline(wks,plot(9),xpts(0:1),ypts(0:1),resl)
  dunck=gsn_add_polyline(wks,plot(10),xpts(0:1),ypts(0:1),resl)
  duncl=gsn_add_polyline(wks,plot(11),xpts(0:1),ypts(0:1),resl)
  ;duncj=gsn_add_polyline(wks,plot(9),xpts(0:1),ypts(0:1),resl)
  ;dunck=gsn_add_polyline(wks,plot(10),xpts(0:1),ypts(0:1),resl)
  ;duncl=gsn_add_polyline(wks,plot(11),xpts(0:1),ypts(0:1),resl)
  ;duncm=gsn_add_polyline(wks,plot(12),xpts(0:1),ypts(0:1),resl)
  ;duncn=gsn_add_polyline(wks,plot(13),xpts(0:1),ypts(0:1),resl)
  ;dunco=gsn_add_polyline(wks,plot(14),xpts(0:1),ypts(0:1),resl)

;;----------------------------------------------------------------------
;; pannel 
;;----------------------------------------------------------------------
  ResP                            = True
;;ResP@gsnMaximize                = True
;;ResP@gsnPaperMargin             = 3 
;;ResP@gsnPanelXWhiteSpacePercent = 1
  ResP@gsnPanelYWhiteSpacePercent = 3 

  gsn_panel (wks,plot, (/4,3/),ResP)


;;----------------------------------------------------------------------
;; output bad site index (0) 
;;----------------------------------------------------------------------
  asciiwrite("kbad.dat",kbad) 

  system("ps2pdf "+ploteps+" "+plotpdf)
  system("mpack -s "+ploteps+" -a "+plotpdf+" k.zhang.iap@gmail.com")

  ;draw(plot)     
  ;frame(wks)


